; +
; NAME:
; gauge_flux_mozart_2D_regrid
;
; PURPOSE:
; Reformats it to mozart grid and writes the file
;
; CATEGORY:
; regridding
;
; CALLING SEQUENCE:
; GAUGE_FLUX_MOZART_2D_REGRID, q_lon, q_lat, q, date, flux_ncname, infilename
;
; INPUTS:
; q_lon = Input longitudes
; q_lat = Input latitudes
; q = Input flux
; date = date from input flux file
; datesec = datesec from input flux file
; flux_ncname = Flux name output from gauge_flux_extract_file_info
; infilename = Input filename
; varinfo = variable information from gauge_flux_extract_file_info
; species = species as string
;
; OUTPUTS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;   Written by:  Ann Stavert 2014
;
; -

Pro GAUGE_FLUX_MOZART_2D_REGRID, q_lon, q_lat, q, date, datesec, flux_ncname, infilename, varinfo, species, TITLE = Title

  ;Create some info for the file
  If keyword_set(TITLE) eq 0 then Title = Species + '_emissions'
  Author='Ann Stavert'
  
  
  ; Calculate the area of each grid cell
  q_area=mr_areagrid(q_lon, q_lat)
  
  ; Need to ensure that the flux variable is in the correct shape
  ; it should be lon x lat x time
  CorrectDims = [ n_elements(q_lon), n_elements(q_lat), n_elements(date)]
  CurrentDims = size(q, /DIMENSIONS)
  ;if N_elements(CurrentDims) eq 2 then CurrentDims = [CurrentDims, 1]
  
  
  if array_equal(CorrectDims, CurrentDims) ne 1 then begin
  
    ; Find the dimensions
    Time_Dim = where(CurrentDims eq n_elements(date), counttime)
    
    Lat_Dim = where(CurrentDims eq n_elements(q_Lat), count)
    Lon_Dim = where(CurrentDims eq n_elements(q_lon), count)
    
    if counttime eq 0 then begin
      ; check if it's only a 2d array i.e. there is no time dimensions
      if size(CurrentDims, /DIMENSIONS) eq 2 then begin
        q_temp = transpose(q, [lon_dim,Lat_dim])
        q_temp = reform(q_temp, [size(q_temp, /DIMENSIONS), 1])
        
        ; set the date variable equal to the central point of the date variable
        date =date[round(N_Elements(date) /2.0)]
        
      endif else begin
        print, 'Can not match array shape to expected array shape'
        return
        
      endelse
      
    endif else begin
      ; time dimension found
      q_temp = transpose(q, [lon_dim,Lat_dim,Time_Dim])
      
    endelse
    q = q_temp
    
  endif
  
  mr_grid_orient, q_lon, q_lat, q, /LONGITUDE
  
  ;change a MOZART-format date ('YYYYMMDD') to Julian date, or visa versa
  ; The date file may be just the month if it's a climatology
  ; need to check for this
  if array_equal(date, indgen(12) + 1) then begin
    ; default to 2000
    print, 'Input date was: [' + strjoin(strtrim(date,2), ', ') + ']'
    print, 'Defaulting to monthly values for the year 2000'
    Date = julday(double(date), 15d, 2000d)
    Date = mr_mozdate(double(Date))
    
  endif
  
  time=mr_mozdate(date)
  time+=double(datesec)/24.d/3600.d
  
  
  ;Get MOZART grid
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  MOZART_fname='/data_rdsf/MOZART/MERRA/MERRA_19x2_20000101.nc'
  lon=mr_ncget('lon', filename=MOZART_fname)
  lat=mr_ncget('lat', filename=MOZART_fname)
  area=mr_areagrid(lon, lat)
  
  
  ;Define output array
  q_regrid=fltarr(n_elements(lon), n_elements(lat), n_elements(time))
  
  ;Martin Shultz's regridding routine.
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  ;Define old grid
  nlon = n_elements(q_lon)
  incr = (max(q_lon) - min(q_lon))/Double(nlon-1)
  ms_lon = Obj_New('MGS_LonVector', $
    start=min(q_lon),incr=incr,nval=nlon, $
    type='regular')
  nlat = n_elements(q_lat)
  incr = (max(q_lat) - min(q_lat))/Double(nlat-1)
  ms_lat = Obj_New('MGS_LatVector', $
    start=min(q_lat),incr=incr,nval=nlat, $
    type='regular')
  o = Obj_New('MGS_RGrid', name='Source Grid', $
    Longitude=ms_lon, Latitude=ms_lat)
  if total(o->getlon() - q_lon) gt 0.1 then stop
  if total(o->getlat() - q_lat) gt 0.1 then stop
  
  ;Define new grid
  LonSize = n_elements(lon)
  lon_incr = (max(lon) - min(lon))/Double(LonSize-1)
  ms_lon = Obj_New('MGS_LonVector', $
    start=min(lon),incr=lon_incr,nval=LonSize, $
    type='regular')
  LatSize = n_elements(lat)
  lat_incr = (max(lat) - min(lat))/Double(LatSize-1)
  ms_lat = Obj_New('MGS_LatVector', $
    start=min(lat),incr=lat_incr,nval=LatSize, $
    type='regular')
  n = Obj_New('MGS_RGrid', name='Target Grid', $
    Longitude=ms_lon, Latitude=ms_lat)
    
  if total(n->getlon() - lon) gt 0.1 then stop
  if total(n->getlat() - lat) gt 0.1 then stop
  
  for ti=0, n_elements(time)-1 do begin
    ;Regrid
    q_regrid[*, *, ti]=regrid(q[*, *, ti], oldgrid=o, newgrid=n, /use_grid, /area_weighted) ;kg/m2/s
    ;      print, string(float(ti)/float(n_elements(time))*100., format='(I)'), '% complete. ' + $
    ;        string(round(total(q_regrid[*, *, ti]*q_area)/total(q[*, *, ti]*area)-1.)*100., format='(I)'), '% error'
    if n_elements(time) gt 1 then $
      print, string(float(ti)/float(n_elements(time))*100., format='(I)'), '% complete. ' + $
      string(round( total(q_regrid[*, *, ti]*area) /total(q[*, *, ti]*q_area) -1.)*100., format='(I)'), '% error' $
      else $
        print, string(float(ti)/float(n_elements(time))*100., format='(I)'), '% complete. ' + $
        string(round( total(q_regrid[*, *]*area)/total(q[*, *]*q_area) - 1.)*100, format='(I)'), '% error'
      
  endfor
  
  ;Clean
  Obj_Destroy, o
  Obj_Destroy, n
  
  ;stop
  
  
  ;CONVERT TO MOLECULES/CM2/S
  ; Extract the units
  AttNames = (varinfo[flux_ncname[0]]).AttName
  Units_Index = where(strmatch(AttNames, 'units', /fold_case))
  Units = ((varinfo[flux_ncname[0]]).Attributes)[Units_Index]
  
  Convert_To_MoleculesCm2Sec, q_regrid, converted, Units, SPECIES =species
  
  
  
  ;Now write in MOZART netCDF file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  ; Write one for each tracer and then one for the sum of the tracers!
  ncName = strjoin((strsplit(InFileName, ".", /EXTRACT))[0:-2],'_') + "_MZT.nc"
  
  fID=ncdf_create(ncName, /clobber)
  
  LonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
  LatDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
  TimeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
  
  ; Create a variable to hold the data
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /DOUBLE)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /DOUBLE)
  dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
  timeID = NCDF_VARDEF(fId,'time', [TimeDID], /DOUBLE)
  qID = NCDF_VARDEF(fId, strcompress(strlowcase(species) + 'emissions', /remove_all), [LonDID, LatDID, TimeDID], /DOUBLE)
  
  ;Create some global attributes
  NCDF_ATTPUT, fid, /Global, "Title", title
  NCDF_ATTPUT, fid, /Global, "Author", author
  NCDF_ATTPUT, fid, /Global, "Created", systime()
  
  ; Create some attributes (about our variables)
  NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
  NCDF_ATTPUT, fId, lonId, "units", "degrees_east"
  
  NCDF_ATTPUT, fId, latId, "long_name", "latitude"
  NCDF_ATTPUT, fId, latId, "units", "degrees_north"
  
  NCDF_ATTPUT, fId, qId, "long_name", Strcompress(strupcase(species) + " emission rate")
  NCDF_ATTPUT, fId, qId, "units", "molecules/cm2/s"
  
  NCDF_ATTPUT, fId, TimeID, "long_name", "Time"
  NCDF_ATTPUT, fId, TimeID, "units", "days since 0000-01-01 00:00:00"
  
  NCDF_ATTPUT, fId, DateID, "long_name", "Date"
  NCDF_ATTPUT, fId, DateID, "units", "YYYYMMDD"
  
  ; Leave definition mode and enter data write mode
  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat
  
  NCDF_VARPUT, fId, dateID, mr_mozdate(time)
  NCDF_VARPUT, fId, timeID, time - julday(1, 121, 1, 0)
  NCDF_VARPUT, fId, qId,  converted
  
  ; Done
  NCDF_CLOSE, fId
  
  print, 'Wrote file : ' + ncName
  
  
END